---
title: "II Article"
author: "Benjamin Toettoe"
date: '2024-09-10'
output:
  html_document:
    df_print: paged
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
```

```{r, include = FALSE}
library(readxl)
library(tidyverse)
library(dplyr)
library(stargazer)
library(lme4)
library(scales)
library(sjPlot)
```

# Importing and cleaning data
```{r, include = FALSE}
dt<-read_excel("C:\\Users\\bento\\OneDrive\\Documents\\UdeM\\Perception Research\\Final Resubmition\\Data.xlsx")
dt$gender <- ifelse(dt$gender == 4, 3, dt$gender)
dt$gender <- ifelse(dt$gender == 3, NA, dt$gender)
dt$pol_lib_con <- ifelse(dt$pol_lib_con == -8, NA, dt$pol_lib_con)
dt$pol_lib_con <- ifelse(dt$pol_lib_con == 99, NA, dt$pol_lib_con)
dt$fbic_bandwidth <- rescale(dt$fbic_bandwidth, to = c(1, 7))
dt$fbic_dependence <- rescale(dt$fbic_dependence, to = c(1, 7))
dt$sat_econ_welbeing<-as.numeric(dt$sat_econ_welbeing)
```

- create means by country
```{r, include = FALSE}
require(dplyr)
head(dt)
dt$country
dt_ag <- dt %>%
  group_by(country) %>%
  summarise(fbic_bandwidth = mean(fbic_bandwidth, na.rm = TRUE),
            fbic_dependence = mean(fbic_dependence, na.rm = TRUE),
            gender = mean(gender, na.rm = TRUE),
            age = mean(age, na.rm = TRUE),
            pol_int = mean(political_interest, na.rm = TRUE),
            pol_lib_con = mean(pol_lib_con, na.rm = TRUE),
            att_country_cn = mean(att_country_cn, na.rm = TRUE),
            att_country_usa = mean(att_country_usa, na.rm = TRUE),
            econ_imp_cn = mean(econ_imp_cn, na.rm = TRUE),
            for_pol_ali_cn = mean(for_pol_ali_cn, na.rm = TRUE),
            educ = mean(educ, na.rm = TRUE),
            sat_econ_welbeing = mean(as.numeric(sat_econ_welbeing), na.rm = TRUE),
            media_soc = mean(media_soc, na.rm = TRUE),
            media_tv = mean(media_tv, na.rm = TRUE),
            media_radio = mean(media_radio, na.rm = TRUE),
            media_newsp = mean(media_newsp, na.rm = TRUE),
            media_online = mean(media_online, na.rm = TRUE),
            econ_pow_cn = mean(econ_pow_cn, na.rm = TRUE),
            covid_us_mil_cn = mean(covid_us_mil_cn, na.rm = TRUE),
            covid_artifical_lab_cn = mean(covid_artifical_lab_cn, na.rm = TRUE),
            political_interest = mean(political_interest, na.rm = TRUE))


```

# Descriptive Statistics (creates Figures 3 and 4)
```{r}
# Calculate the average level of perceived economic importance for each country
average_per_country <- aggregate(dt$econ_imp_cn, by = list(dt$country), FUN = function(x) mean(x, na.rm = TRUE))
colnames(average_per_country) <- c("country", "average_econ_imp_cn")

# Sort the data frame by the average value
average_per_country <- average_per_country[order(average_per_country$average_econ_imp_cn), ]

# Calculate standard errors for the error bars
standard_errors <- tapply(dt$econ_imp_cn, dt$country, FUN = function(x) sd(x, na.rm = TRUE) / sqrt(sum(!is.na(x))))
standard_errors_df <- data.frame(country = names(standard_errors), standard_error = unlist(standard_errors))

# Combine the data frames for the average values and standard errors
average_per_country <- merge(average_per_country, standard_errors_df, by = "country")

plot_1 <- ggplot(average_per_country, aes(x = average_econ_imp_cn, y = reorder(country, -average_econ_imp_cn))) +
  geom_point(stat = "identity", fill = "lightblue") +
  geom_errorbarh(aes(xmin = average_econ_imp_cn - standard_error, xmax = average_econ_imp_cn + standard_error), height = 0.2) +
  labs(y = "Country", x = "Perceived Economic Importance of China", title = "") +
  theme(axis.text.y = element_text(hjust = 1)) + theme_bw()


### now let's do fbic_dependence
# Calculate the average level of fbic_dependence for each country
average_fbic_dependence <- aggregate(dt$fbic_dependence, by = list(dt$country), FUN = mean)
colnames(average_fbic_dependence) <- c("country", "average_fbic_dependence")

# Sort the data frame by the average value
average_fbic_dependence <- average_fbic_dependence[order(average_fbic_dependence$average_fbic_dependence), ]

# Calculate standard errors for the error bars
standard_errors_fbic <- tapply(dt$fbic_dependence, dt$country, FUN = function(x) sd(x) / sqrt(length(x)))
standard_errors_fbic_df <- data.frame(country = names(standard_errors_fbic), standard_error_fbic_dependence = unlist(standard_errors_fbic))

# Combine the data frames for the average values and standard errors
average_fbic_dependence <- merge(average_fbic_dependence, standard_errors_fbic_df, by = "country")

# Create the graph with error bars (points) for fbic_dependence
plot_2 <- ggplot(average_fbic_dependence, aes(x = average_fbic_dependence, y = reorder(country, -average_fbic_dependence))) +
  geom_bar(stat = "identity") +
  geom_errorbarh(aes(xmin = average_fbic_dependence - standard_error_fbic_dependence, xmax = average_fbic_dependence + standard_error_fbic_dependence), height = 0.2) +
  labs(y = "Country", x = "Economic Dependence on China (FBIC)", title = "") +
  theme(axis.text.y = element_text(hjust = 1)) + theme_bw()

plot_1
plot_2

```

plot the two variables together into one scatterplot (creates Figure 5)
```{r}
library(lme4)
library(sjPlot)
library(ggplot2)


# Calculate the average level of perceived economic importance for each country
average_econ_imp_cn <- aggregate(dt$econ_imp_cn, by = list(dt$country), FUN = function(x) mean(x, na.rm = TRUE))
colnames(average_econ_imp_cn) <- c("country", "average_econ_imp_cn")

# Calculate the average level of fbic_dependence for each country
average_fbic_dependence <- aggregate(dt$fbic_dependence, by = list(dt$country), FUN = mean)
colnames(average_fbic_dependence) <- c("country", "average_fbic_dependence")

# Combine the two data frames
combined_data <- merge(average_econ_imp_cn, average_fbic_dependence, by = "country")

# Create the scatter plot
scatter_plot<-ggplot(combined_data, aes(x = average_fbic_dependence, y = average_econ_imp_cn)) +
  geom_point(color = "grey", size = 3) +
  geom_smooth(method = 'lm')+
  geom_text(aes(label = country), hjust =- 0.25, vjust = -0.1) +
  labs(x = "Perceived Economic Importance of China", 
       y = "Economic Dependence on China (FBIC)",
       title = "") +
  theme_minimal()

# Print the scatter plot
print(scatter_plot)

```

# Analysis
- Create aggregate models (Creates Figure 7 and Appendix 5)
```{r, echo = FALSE}
library(car)

a1 <- lm(for_pol_ali_cn ~ econ_imp_cn + fbic_dependence, data = dt_ag)
a2 <- lm(for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + age + gender, data = dt_ag)
a3 <- lm(for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + pol_int + pol_lib_con, data = dt_ag)
a4 <- lm(for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + age + gender + pol_lib_con + pol_int,  data = dt_ag)

stargazer(list (a1, a2, a3, a4), title = "Aggregate Regression Output",
          star.cutoffs = c(.05, .01, .001),
          align = TRUE, type = 'text',
          omit.table.layout = "n")


correlation_coefficient <- cor(dt_ag$econ_imp_cn, dt_ag$fbic_dependence, use = "complete.obs")
print(correlation_coefficient)


plot_model(a4, show.values = TRUE, value.label = format(p_values, digits = 3),
           theme = sjPlot::theme_sjplot(base_size = 40))
dev.off()


```

- Multi-level models (Creates Figure 6 and Appendix 4)
```{r, echo = FALSE}
dt$country<-as.factor(dt$country)
dt$country <- relevel(dt$country, ref = 'Colombia') # baseline country is colombia

m1 <- for_pol_ali_cn ~ econ_imp_cn + fbic_dependence +  (1 + econ_imp_cn| country)
m2 <- for_pol_ali_cn ~ econ_imp_cn + fbic_dependence +  (1| country) # only fixed effect
m3 <- for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + age + gender + (1 | country)
m4 <- for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + pol_lib_con + political_interest  + (1 | country)
m5 <- for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + age + gender  + pol_lib_con + political_interest + (1 | country)


mreg1 <-lmer(m1, data = dt)
mreg2<-lmer(m2,data = dt)
mreg3<-lmer(m3,data= dt)
mreg4<-lmer(m4,data= dt)
mreg5<-lmer(m5,data = dt)


stargazer(list (mreg1, mreg2, mreg3, mreg4, mreg5), title = "Regression Output",
          star.cutoffs = c(.05, .01, .001),
          align = TRUE, type = 'text',
          omit.table.layout = "n")

vif(mreg5)

correlation_coefficient <- cor(dt$econ_imp_cn, dt$fbic_dependence, use = "complete.obs")
print(correlation_coefficient)

plot_model(mreg5, show.values = TRUE, value.label = format(p_values, digits = 3),
           theme = sjPlot::theme_sjplot(base_size = 40))
dev.off()

```

- Naive Mediation Analysis (Creates Figure 8 and first table in Appendix 6)
```{r}
set.seed(123456)
library(multilevel)
library(mediation)
fitM<-lmer(econ_imp_cn ~ fbic_dependence +  (1| country),data=dt)
fitY<-lmer(for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + (1 | country),data = dt)


#now we are doing the actual mediation analysis
within_effects <- fixef(fitM)
between_effects<-ranef(fitY)
between_effects<-between_effects$country[,'(Intercept)']
med_result <- mediation::mediate(fitM,fitY, treat = "fbic_dependence", mediator = "econ_imp_cn")
summary(med_result)
plot(med_result) 
modelsummary:: modelsummary(fitY, star = T)

summary(fitY)
stargazer(fitY,star.cutoffs = c(.05, .01, .001), type = 'text',report=("vc*p"))

#ACME: average causal mediaiton effects (indirect effect of IC on DV)
#ADE: average direct effect

library(diagram)

dt_md<- c(0, "'.593*'",0,
          0,0,0,
          "'.673***'","'.291***'",0)
M<- matrix(nrow = 3, ncol = 3, byrow = T, data = dt_md)
plot<-plotmat(M, pos = c(1,2),
              name = c('Perceived Chinese 
Economic Importance', 'Actual Dependence 
on China','Alignment with China'),
              box.type = 'rect', box.size = 0.15, box.prop = 0.25, curve = 0)
```


- Mediation Analysis with Controls (Creates second table in Appendix 6)
```{r, eval = FALSE}
set.seed(123456)
library(multilevel)
library(mediation)
fitMM<-lmer(econ_imp_cn ~ fbic_dependence + age + gender + political_interest + pol_lib_con + (1| country), data=dt)
summary(fitMM)
fitYY<-lmer(for_pol_ali_cn ~ fbic_dependence + econ_imp_cn +  age + gender + political_interest + pol_lib_con + (1 | country),data = dt)

within_effects_2 <- fixef(fitMM)
between_effects_2<-ranef(fitYY)
between_effects_2<-between_effects_2$country[,'(Intercept)']
med_result_2 <- mediation::mediate(fitMM,fitYY, treat = "fbic_dependence", mediator = "econ_imp_cn")
summary(med_result_2)
plot(med_result_2) 
modelsummary:: modelsummary(fitYY, star = T)

#ACME: average causal mediaiton effects (indirect effect of IC on DV)
#ADE: average direct effect

library(diagram)
dt_md<- c(0, "'.101***'",0,
          0,0,0,
          "'.0.127***'","'.252***",0)

dt_md<- c(0, "'.101***'",0,
          0,0,0,
          "'.187***'","'.187**'",0)
M<- matrix(nrow = 3, ncol = 3, byrow = T, data = dt_md)
plot<-plotmat(M, pos = c(1,2),
              name = c('Perceived Chinese 
Economic Importance', 'Actual Dependence 
on China','Alignment with China'),
              box.type = 'rect', box.size = 0.15, box.prop = 0.25, curve = 0)
```


- Multi-level model with individualistic interaction terms (Creates Appendix 7)
```{r, echo = FALSE}

i1 <- for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + age + gender + pol_lib_con + political_interest + educ + sat_econ_welbeing + (1 | country)
i2 <- for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + age + gender + pol_lib_con + political_interest + educ + sat_econ_welbeing + fbic_dependence:educ + (1 | country)
i3 <- for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + age + gender + pol_lib_con + political_interest + educ + sat_econ_welbeing + fbic_dependence:sat_econ_welbeing +(1 | country)
i4 <- for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + age + gender + pol_lib_con + political_interest + educ + sat_econ_welbeing + econ_imp_cn:educ + (1 | country)
i5 <- for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + age + gender + pol_lib_con + political_interest + educ + sat_econ_welbeing + econ_imp_cn:sat_econ_welbeing + (1 | country)


ireg1 <-lmer(i1, data = dt)
ireg2<-lmer(i2,data = dt)
ireg3<-lmer(i3,data= dt)
ireg4<-lmer(i4,data= dt)
ireg5<-lmer(i5,data = dt)


stargazer(list (ireg1, ireg2, ireg3, ireg4, ireg5), title = "Regression Output",
          star.cutoffs = c(.05, .01, .001),
          align = TRUE, type = 'text',
          omit.table.layout = "n")



plot_model(ireg4, show.values = TRUE, value.label = format(p_values, digits = 3),
           theme = sjPlot::theme_sjplot(base_size = 40))
dev.off()


plot_model(ireg5, show.values = TRUE, value.label = format(p_values, digits = 3),
           theme = sjPlot::theme_sjplot(base_size = 40))
dev.off()


```

marginal effect plots (Creates Figure 9 and 10)
```{r}
custom_labels <- c("Very Dissatisfied", "Neutral", "Very Satisfied")
plot_model(ireg2,type = 'pred',terms = c("fbic_dependence","educ"))
plot_model(ireg3,type = 'pred',terms = c("fbic_dependence","sat_econ_welbeing[1,4,7]"),
           legend.title = "Economic Wellbeing",
           axis.title = c("FBIC Dependence", "Predicted Value of Foreign Policy alignment with China"),
           labels = custom_labels)+labs(title = "")+scale_color_manual(values = c("red", "blue", "green"),labels = custom_labels)

##
i4 <- for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + age + gender + pol_lib_con + political_interest + educ + sat_econ_welbeing + econ_imp_cn:educ  + (1 | country)
ireg4<- lmer(i4,data=dt)
stargazer(ireg4,type = 'text')

i5 <- for_pol_ali_cn ~ econ_imp_cn + fbic_dependence + age + gender + pol_lib_con + political_interest + educ + sat_econ_welbeing + econ_imp_cn:sat_econ_welbeing + (1 | country)
ireg5<- lmer(i5,data = dt)
plot_model(ireg4,type = 'pred',terms = c("econ_imp_cn","educ"))
custom_labels <- c("Very Dissatisfied", "Neutral", "Very Satisfied")
plot_model(ireg5,type = 'pred',
           terms = c("econ_imp_cn","sat_econ_welbeing[1,4,7]"),
           legend.title = "Economic Wellbeing",
           axis.title = c("Econ Imp CN", "Predicted Value of Foreign Policy alignment with China"))+labs(title = "")+scale_color_manual(values = c("red", "blue", "green"),labels = custom_labels)

```




Disconnect Analysis (Creates Appendix 8 and 9)
```{r, echo = FALSE}
library(dplyr)
dt$disconnect <- dt$econ_imp_cn - dt$fbic_dependence
dt$disinform <- (dt$covid_us_mil_cn + dt$covid_artifical_lab_cn)/2

overestimate <- dplyr::filter(dt, disconnect > 0)
underestimate <- dplyr::filter(dt, disconnect < 0)

d1 <- lm(disconnect ~ disinform, data = overestimate)
d2 <- lm(disconnect ~ disinform + age + gender, data = overestimate)
d3 <- lm(disconnect ~ disinform + pol_lib_con + political_interest, data = overestimate)
d4 <- lm(disconnect ~ disinform + age + gender + pol_lib_con + political_interest, data = overestimate)

d5 <- lm(disconnect ~ disinform, data = underestimate)
d6 <- lm(disconnect ~ disinform + age + gender, data = underestimate)
d7 <- lm(disconnect ~ disinform + pol_lib_con + political_interest, data = underestimate)
d8 <- lm(disconnect ~ disinform + age + gender + pol_lib_con + political_interest, data = underestimate)

stargazer(list (d1, d2, d3, d4), title = "Regression Output for those overestimating China's importance",
          star.cutoffs = c(.05, .01, .001),
          align = TRUE, type = 'text',
          omit.table.layout = "n")

stargazer(list (d5, d6, d7, d8), title = "Regression Output for those Underestimating China's importance",
          star.cutoffs = c(.05, .01, .001),
          align = TRUE, type = 'text',
          omit.table.layout = "n")

d9 <- lm(disconnect ~ media_soc + media_newsp + media_tv + media_radio + media_online, data = overestimate)
d10 <- lm(disconnect ~ media_soc + media_newsp + media_tv + media_radio + media_online + age + gender, data = overestimate)
d11 <- lm(disconnect ~ media_soc + media_newsp + media_tv + media_radio + media_online + pol_lib_con + political_interest, data = overestimate)
d12 <- lm(disconnect ~ media_soc + media_newsp + media_tv + media_radio + media_online + age + gender + pol_lib_con + political_interest, data = overestimate)

d13 <- lm(disconnect ~ media_soc + media_newsp + media_tv + media_radio + media_online, data = underestimate)
d14 <- lm(disconnect ~ media_soc + media_newsp + media_tv + media_radio + media_online + age + gender, data = underestimate)
d15 <- lm(disconnect ~ media_soc + media_newsp + media_tv + media_radio + media_online + pol_lib_con + political_interest, data = underestimate)
d16 <- lm(disconnect ~ media_soc + media_newsp + media_tv + media_radio + media_online + age + gender + pol_lib_con + political_interest, data = underestimate)

stargazer(list (d9, d10, d11, d12), title = "Regression Output for those overestimating China's importance",
          star.cutoffs = c(.05, .01, .001),
          align = TRUE, type = 'text',
          omit.table.layout = "n")

stargazer(list (d13, d14, d15, d16), title = "Regression Output for those Underestimating China's importance",
          star.cutoffs = c(.05, .01, .001),
          align = TRUE, type = 'text',
          omit.table.layout = "n")

library(sjPlot)
modelsummary::modelsummary(d12, star = T)
plot_model(d12, show.values = T,  value.label = format(p_values, digits = 3))

library(sjPlot)
modelsummary::modelsummary(d16, star = T)
plot_model(d16, show.values = T,  value.label = format(p_values, digits = 3))
```